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We propose a spatially continuous and temporally dis- 
crete model for pattern formation in vertically vibrated gran- 
ular layers. The grain transfer and the grain mobility transi- 
tions are introduced qualitatively, but explicitly. This model 
reproduces experimentally observed f/2 patterns, includ- 
ing squares in high-mobility conditions and stripes in low- 
mobility conditions. Hexagons are obtained when a high- 
mobility condition and a low-mobility condition are alter- 
nately applied. We discuss this hexagon formation in connec- 
tion with the formation of squares and stripes. In addition 
to the experimentally observed localized oscillon structures, 
this model also exhibits localized hexagonal structures. These 
localized hexagonal structures grow by nucleating oscillons of 
equal phase. 

PACS numbers: 45.70.Qj, 05.45.-a, 83.10.Ji, 45.70.-n 



I. INTRODUCTION 

Spatiotemporal patterns are found in numerous phys- 
ical, chemical, and biological systems, and have drawn 
much interests in recent years [^]j2]. Like fluid sys- 
tems (Faraday systems) , familiar but poorly under- 
stood granular systems ||^ show various patterns includ- 
ing stripes, squares, hexagons and kinks when vertically 
vibrated Novel two-dimensional localized struc- 

tures, termed 'oscillons', were also observed in these sys- 
tems 0,^. Since the pattern formation phenomena in 
vertically vibrated granular systems were first reported, 
many models have been proposed to explain the mecha- 
nism pO|-pO| . Molecular dynamics simulations quan- 
titatively reproduce every observed patterns except os- 
cillons and kinks, but do not give the essentials of the 
mechanism. Phenomenological models ||Tl|-p^ explain 
some of the phenomena by considering the grain trans- 
fer with a few qualitative features. Especially, the model 
in Ref. jlj] reproduces stripes, squares, and hexagons 
by using two simple features, i.e., randomizing impacts 
and inelastic collisions. The continuum coupled map 
model yields a phase diagram similar to that of 

the experiments with the particular choice of the map 
and the introduction of a second length scale and ex- 
plains the universal bifurcation sequence related to the 
interaction between temporal period doubling and spatial 
pattern formation. The results obtained by using the or- 
der parameter equation indicate that the localized 
structures are not specific to granular systems p6|,p7|. 



Although these models give the key to the mechanism, 
hexagon formation has not been properly treated. The 
models of Refs. 
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, |18| | exhibit hexagons, but do not 
elucidate the hexagon forming conditions. Other models 
reproduce only squares Q or only stripes besides 
oscillons. 

In this paper, by incorporating the grain transfer 
p2|~p^ qualitatively but explicitly, we propose a model 
that can cover all patterns including squares, stripes, 
hexagons, and oscillons. Hexagon formation will be dis- 
cussed in connection with squares and stripes. 

Previously, we proposed a simple model which 
gave some insights into the dynamics of stripes, oscillons, 
and kinks. The following two features were assumed to 
study stripes and oscillons |p^ : 

(i) A local excitation organizes only through interac- 
tion with excitations within some distance; 

(ii) The effective interaction or flow is approximately 
determined by the difference between the excitations. 
Additionally, a period-doubling bifurcation nature of 
the local dynamics was needed to investigate the kinks 
|]6|-^, ^. With these features and a Poincare map-like 
idea ]21[ | , we simplified the dynamics of the layer and 
constructed a model which maps the pattern at the time 
of the layer-plate collision to the pattern at the time of 
the next collision. 

However, room exists for more generalizations of the 
model, for it does not show squares or hexagons yet. In 
this regard, some authors P, pl] , p^p^ , |2^ argue that mass 
conservation has an important role in pattern formation 
and grain mobility contributes to the selection of squares 
and stripes. In experiments |^,^, squares arise at low 
forcing frequencies, and at those frequencies, the hori- 
zontal mobility is large due to the large dilation of the 
layer and the small dissipation. However, stripes are ob- 
served at high frequencies, where the mobility is small 
due to the small dilation and the large dissipation. 

In Sec. II, we propose a modified model by introducing 
mass conservation and grain mobility transitions. The re- 
sults from numerical simulations are presented and dis- 
cussed in Sec. HI. 



II. MODEL 

In this modified model, we use densities instead of lo- 
cal heights to consider the grain transfer explicitly. The 
density (t„ (r ) denotes how many grains are piled up on 
a position r in a two-dimensional space at the nth time 
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step. As mentioned above, the nth time step corresponds 
to the nth coUision between the layer and the plate. Note 
that in the case of //2 patterns, the densities behave as 
the local heights of the patterns we observe. 

This model maps the densities at the nth time step to 
those at the (n + l)th time step through following two 
stages: 



The first stage : (r'n^i^f ) = (T„(r' ) + Aa{f ), 



(1) 



The second stage : dt<T{r,t) — DV ±^(7{f,t), (2) 
or equivalently, 

) = e-^*'''=V'„+i(^), (3) 



where Aa{f ) is the increment or the decrement in the 
density at r, D is a diffusion coefficient, td is the effec- 
tive time duration of diffusion, and y±^ = dl + dy. In 

Eq. (H), (T„(fc ) denotes the spatial Fourier transform of 

The two stages are similar to those of the model pro- 
posed in Ref. [ p^ , but here we can control the grain mo- 
bility. Also, as our original model [ psljl^ ], this model 
shares a common spirit with the model of Ref. p8| , p^ : 
these models are spatially continuous, but temporally 
discrete. 

The First Stage. Since there is a limit on the dis- 
tance a grain can travel during the time interval be- 
tween the layer-plate collisions, we take the inside of a 
circle with radius R centered at r in two-dimensional 
space as an interaction region. Some of the grains 
at this position f will flow into the interaction region 
M{f ) = {r"' I < \r' -f\<R}. The radius R can 
be determined by the horizontal velocity, the free flight- 
time, and the collisions with other grains. For simplicity, 
we set R to the same constant for all f's. 

The number of grains transferred from one position 
to another depends largely on the lateral velocity of the 
grains and the free-flight time during which the lateral 
transfer of grains occurs 13|. The lateral velocity is 
approximately proportional to both the layer-plate colli- 
sion velocity and the slope of the free surface of the layer 
Ip^ ; the free flight-time is roughly proportional to 1//. 

From these facts, we consider the quantities 
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^G(.,,0-...(.-))M,(f;^), (4, 



S„(r)=/ dr' Snir^r'), r' eM<{r), (5) 

iA/'<(r) 

where A/'<(r ) = {r' e Af{f )\an{r') < cr„(f )}. The 
quantity 5„(r f') in Eq. (0) is the maximum num- 
ber of grains that can be transferred from r to t^' where 
the density is lower than f if there are enough grains at 
position r. We call this maximum number the demand 



from f' to f. The actual number transferred is limited 
by the density at f. For this purpose, S]„(r ), the total 
demand to f, is given in Eq. (^. In this model, the pa- 
rameter a acts like the layer-plate collision velocity (or 
the relative collision acceleration) which determines the 
onset of patterns [|6|jl^. The dependence of the trans- 
ferred number on the slope of the free surface is phe- 
nomenologically represented by the function G, which 
depends only on the difference between the densities at 
the two positions. Due to the avalanche-type flow p|,p5|, 
the number of grains transferred is naturally expected to 
increase faster than linearly as the difference increases, 
so we choose a slope- increasing function for G [G" > 0; 
see Fig. 1 (a)]. As in our previous model |T^,|l^, this 
slope-increasing property of the function G makes the 
hysteresis more evident and strong: As patterns grow 
from the small fluctuations in the densities, the density 
differences increase more, and the effect of this coupling 
function becomes more important. Therefore, this effect 
makes the quantities in Eqs. (Q) and (||) large enough 
to support the excitations even if a is lowered below the 
critical values where the transitions from flat to patterns 
occur (see Fig. 4 to be discussed). The function 
controlling the mobility of the grains at a position r is a 
function of the status of the position r in the interaction 
region Af{r ). e is a control parameter. The status of a 
position r is given by the quantity S'(A/'<(r ))/S'(A/'(r )), 
where S'(A/<(r )) is the area with densities lower than 
(T„(r ) in the interaction region Af{r), and S'(A/'(r)) is 
the total area of the interaction region. Thus, positions 
with relatively higher densities have higher status. For 
example, a peak in the layer corresponds to the highest 
status and a valley to the lowest one. Experimentally, 
there are qualitative changes in the overall mobility as 
the forcing frequency changes. The mobility is high at 
low frequencies while at high frequencies the grains are 
less mobile. Without M^, in this model, the grains at po- 
sitions with high status are more apt to move than those 
at positions with low status due to the large values of the 
quantities Sn{f —> r*') and S„(r ). We can control the 
relative mobilities between these positions by giving dif- 
ferent values to them according to their status. If we 
give high values to positions with low status [see Fig. 
I (b), e < 0], the mobilities at these positions increase. 
In such a case, most of the grains are apt to move. Con- 
trary to this case, if we give low values to positions 
with low status [see Fig. I (b), e > 0], the grains at these 
positions are even more immobile. The parameter e in 
our model roughly plays a role analogous to that of the 
forcing frequency in the experiments. The denominator 
S{Af{r)) inEq. (|) is introduced to diminish the a value 
dependence on the radius of the interaction region J\f{r). 

Since there are limited number of grains at the position 
r and we want to conserve the mass, satisfying equally 
the demands from f"' , we distribute the limited number 
of grains to A/< (r ) in two possible ways: 

Case 1: If the total demand S„(r ) is less than the 
density at r, the outflow from ftof'e A/<(f ) is set to 



2 



the demand 5„(r ^ f'). 

Case 2: If the total demand S„(r ) is greater than 
the density at r, the total outflow from r is set to the 
density at f, and the outflow is distributed to A/< (r ) in 
proportion to the demand (5n(r — > r'). 
Note that grains are more actively transferred in Case 2 
than in Case 1, for in such a case, all the grains at those 
positions leave their positions. 

In this manner, we can calculate the inflow and the 
outflow of the grains for each position at this first stage. 

The Second Stage. In the second stage, to include 

the diffusion effect due to the collisions between grains 
and between the layer and the plate [ p"3yi^ ] , we apply Eq. 

which automatically conserves mass. This process re- 
laxes the high-density gradient occurring within the size 
of the characteristic length scale Aq = 2R. Since the net 
diffusion effect is determined by the parameter Dtd, we 
choose Dtd ~ (x^)~^ such that the Fourier components 
of the density distribution with wave vectors greater than 
fco = 1^ are sufficiently suppressed. 



III. NUMERICAL RESULTS AND DISCUSSION 

For the numerical simulations, we choose the functions 
(see Fig. 1) 



G{x) — X -\- cx^ , c > 0, X > 0, 
M^{x) = < a: < 1. 



(6) 
(7) 



These functions are the simplest ones with the required 
properties. G{x) in Eq. is a monotonically increasing 
function with a slope-increasing property (G" > 0). 
in Eq. is monotonically increasing for e > 0, but 
monotonically decreasing for e < 0. Other functions with 
the same properties do not give qualitatively different 
results. We fix c = 0.7 in the following simulations for the 
stability of oscillon structures. Higher (lower) c values 
give larger (smaller) hysteresis. 

We performed numerical simulations on a 128 x 128 
grid with spatial extension L — 128 and periodic bound- 
ary conditions. Initially, the densities are distributed 
around the mean value ctq = 10 with small uniform ran- 
dom deviations within [— lO^'^, lO^'^]. 

Square and stripe patterns obtained using this model 
are shown in Figs. 2 and 3. These patterns oscillate 
with //2 as in experiments [^0]; that is, they return 
to the same configuration after 2 periods of the forcing 
(in this case, one time step corresponds to one period). 
Figure 4 shows the stability region for each pattern as 
a function of a and e. There are two transition regions 
between squares and stripes. At lower a's, the transi- 
tion from squares to stripes occurs as the parameter e 
increases from negative to positive values, and for posi- 
tive e's, the transition from stripes to squares occurs as a 
increases. Both of these transitions correspond to quali- 
tative changes in the grain mobility. As expected, squares 



arise for e < or sufficiently large a, where the quanti- 
ties in Eqs. (||) and (||) are large enough that grains are 
actively transferred: All the grains at positions of more 
than 85% of the entire space leave their positions accord- 
ing to Case 2 at the first stage of our model. On the 
contrary, when stripes appear, the percentage decreases 
below 80%. For the same reason, greater c values for 
G{x) in Eq. (^ lowers the transition points from stripes 
to squares for positive e values. 

Figure 3 also reflects the grain mobility transitions. 
As pointed out in Ref . ||l^ , the densities are sharply clus- 
tered in the square patterns (the high-mobility situation) 
and more diffuse in the stripe patterns (the low-mobility 
situation); the peaks in the square pattern are higher 
than those in the stripe patterns (see Fig. 3). 

In the hysteretic region, this model also exhibits 
period-2 oscillons and the oscillon structures such as 
dipolcs, triangular tetramers, and chains (not shown), 
which are observed experimentally [^j7|. When a is in- 
creased above the upper stability limit, the structures 
grow by nucleation of oscillons (e < in this model, 
low / in experiments) or by elongation of oscillons and 
nucleation of waves (e > in this model, / > 30 in ex- 
periments) [[tIPIP^. 

Now, we consider the hexagon formation. Hexagons 
are also period-2 patterns like squares and stripes, but 
there is a difference. Two distinct phases - a set of iso- 
lated peaks on a triangular lattice and a hexagonal cell 
- appear in alternate cycles. Experimentally, hexagons 
arise when effective two-frequency forcing is applied in- 
ternally or externally In this situation, the layer 
has two different free-flight times and, consequently, col- 
lides with the plate alternately with two different col- 
liding conditions. In the absence of any clear identifi- 
cation of the colliding conditions |^,|l9|, from the two 
phases of hexagons, we conjecture that one is a high- 
mobility condition, and the other is a low-mobility con- 
dition. To determine the validity of our conjecture, we 
apply a high-mobility condition (a square-forming con- 
dition) in one step and in the next step a low-mobility 
condition (a stripe-forming condition) with the same ra- 
dius. This yields hexagonal patterns as shown in Figs. 5 
and 6. 

As in experimental results ^, hexagonal patterns 
obtained in this model show two different phases - a peak 
phase and a cellular phase (Fig. 5). Note that the peak 
phase has localized density peaks like squares and always 
appears after applying the high-mobility condition [see 
Fig. 3(a), Fig. 6(a), Fig. 5(a)]. Similarly, note also that 
the cellular phase has ridges like stripes and always arises 
after applying the low- mobility condition [see Fig. 3(b), 
Fig. 6(b), Fig. 5(b)]. From these facts, we can regard the 
hexagon formation as a compromise between squares and 
stripes. Our results are consistent with the occurrence of 
hexagons in the transition region between squares and 
stripes in Fig. 2 of Ref. [|l^. However, the hexagon 
formation described by the order parameter equation p^ ] 
is different from ours. In our model, grain transfer is 
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crucial for hexagon formation, but in the case of Ref. 
grain transfer is treated irrelevant to hexagon formation. 

Since the hexagons in this model are formed directly 
from the flat state and since the transitions are hysteretic, 
we have parameter regions with coexistence of hexagons 
and flat. In that region, we obtain stable localized hexag- 
onal structures as shown in Fig. 7. These structures 
may be regarded as oscillon structures. Actually, if the 
parameter a is raised above the upper stability limit, 
these hexagonal structures grow into a global hexago- 
nal pattern by nucleating oscillons of equal phase at the 
hexagonal lattice sites (see Fig. 8) as the oscillon struc- 
tures at low forcing frequency grow through square-like 
structures into squares and at high frequency through 
stripe-like structures into stripes Although such 

localized structures have not been reported in granular 
systems, similar structures have been found in other sys- 
Especially, the localized hexagonal oscil- 
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tems 

Ion structures found in the two-frequency forcing exper- 
iments with Newtonian fluids seem to be closely re- 
lated to the structures found in our model. The growing 
phenomena reported in Refs. |^,^,^ are also similar 
to the one described above. In contrast, the hexagonal 
patch in Ref. p8) grows by a six-faceted, moving front, 
and the three independent roll states in superposition 
give the hexagonal structure. 



IV. CONCLUSION 

In summary, considering the grain transfer qualita- 
tively but explicitly, we have constructed a model for 
pattern formation in vertically vibrated granular layers. 
Our model maps the density distribution at the time of 
the layer-plate collision to that at the time of the next col- 
lision in a continuous two-dimensional space. This model 
exhibits squares in high-mobility conditions and stripes 
in low-mobility conditions. To elucidate the hexagon 
forming situation, we simulated period doubling situation 
by applying a high-mobility condition (a square-forming 
condition) and a low-mobility condition (a stripe-forming 
condition) alternately and observed hexagons. Our sim- 
ulations show that the peak phase of hexagons appears 
after the high-mobility condition and the cellular phase 
after the low-mobility condition. The peak phase has lo- 
calized density peaks like squares, and the cellular phase 
has ridges like stripes. From these facts, we conclude that 
in vertically vibrated granular layers hexagon formation 
is a compromise between squares and stripes. Finally, the 
model presented here also reproduces the experimentally 
observed oscillon structures and their dynamics. In ad- 
dition, we observed localized hexagonal structures and 
their growth by nucleation. Our model may be use- 
ful for the studies of multifrequency forcing experiments 
26|j27|j3^ and other pattern forming systems PlJSS 
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Figure Captions 

Fig. 1. Functions in Eq. (^. (a) G{x) is a function of 
the density difference between two positions, (b) 
Mg is a function of the status (see text for defini- 
tion), the parameter e controls the overall grain 
mobility. 

Fig. 2. Shadow graph of the patterns obtained numer- 
ically in our model. (White denotes high density.) 
(a) Period-2 square pattern at a = 3, e = —3, R = 
6, and Dtd — 3.65. (b) Period-2 stripe pattern at 
a = 11, e = -1-3, R — 6, and Dtd — 3.65. 

Fig. 3. Three-dimensional perspective of the densities 
in (a) the square and (b) the stripe patterns of Fig. 
2. 

Fig. 4. Phase diagram of our model as a function of 
a and e with R = 6 and Dtd = 3.65. Closed 
(open) squares denote transitions from flats (pat- 
terns); closed triangles and circles indicate the tran- 
sition points. The numbers near symbols are the 
corresponding a values. 

Fig. 5. Hexagonal pattern obtained by alternately ap- 
plying a high-mobility condition and a low-mobility 
condition, (a) Peak phase {n — odd) after apply- 
ing the high-mobility condition (ai = 5,ei = —3) 
and (b) cellular phase (n — even) after apply- 
ing the low-mobility condition {a2 — 5,e2 = +3). 
i?i = i?2 - 6 and {Dtd)i = {Dtd)2 - 3.65. 

Fig. 6. Three-dimensional perspective of the densities 
in the hexagonal pattern of Fig. 5 in (a) the peak 
phase and (b) the cellular phase. 

Fig. 7. Localized hexagonal structure in the hysteresis 
region obtained by applying the same conditions as 
in Fig. 5, except for ai = a2 = 1.1. 

Fig. 8. Growth of hexagonal structure by nuclcation. 
The a values are raised from 1.1 in Fig. 7 to 1.5. 
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